arXiv:l506.05319vl [math.ST] 17 Jun 2015 


Cumulants of products of Normally distributed 
random variables 

Clarence KALITSI and Jan VRBIK 
Department of Mathematics 
Brock University, 500 GLenridge Ave. 

St. Catharines, Ontario, Canada, L2S 3Al 


Abstract 

To find moments of various estimators related to Autoregressive mod¬ 
els of Statistics, one first needs the cumulants of products of two Normally 
distributed random variables. The purpose of this article is to derive the 
corresponding formulas, and extend them to products of three or more 
such variables. 


1 Introduction 

The formulas presented in this article are crucial for finding moments (and, sub¬ 
sequently, approximate distributions) of various parameter estimators related to 
AR(k) models (see [1] and El). 


1.1 Multivariate Normal distribution 


Assume that Xi, V 2 ,... are centralized (having zero mean), Normally distributed 
random variables. Their moment-generating function is given by 


exp 


(^) 
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where V is the corresponding variance/covariance matrix. Based on m , we 
can easily find the expected value of a product of any number of such random 
variables (which defines the corresponding moment), getting zero when this 
number is odd, and 


E {XiX2...X2k) = A{1,2,. 
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when this number is even. The summation is over all possible (2 ^ 2 )/^- w^iys 
of PAIRING the n indices, regardless of the order (both within and between) 
of the resulting pairs. Thus, for example, {1,2}, {3,4} and {4,3}, {2,1} are 
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considered identical and appear only once, but {1,2}, {3,4} and {1,3}, (2,4} 
are distinct. 

Notationally (and rather symbolically), this has been indicated by 

ltV2k 

where i represents the {ii,i 2 }, {isjU},--- indices and V 2 k the set of all such 
selections. 

Example 

A*{1.2,3,4} = {XiX2X3Xi) = Vi_2V3^4 + Vi,3V2,4 + hi,4V2,3 til (3) 
Note that the resulting formulas are fully general; m can be read as 

^{ilb2 J 3 J 4 } ~ y?! J2^3J4 + y?! J3^2 J4 + ^'104^2 J3 

but using specific integers simplifies the notation; also, the formulas allow any 
duplication of indices, e.g. 

M{ 14 , 2 , 2 } = td,lV2,2 + 2Vi2 

It is important to realize that, for random variables with zero mean (the vari¬ 
ables of this section), there is no difference between simple and central moments; 
our fi thus stands for either. 


1.2 Joint cumulants 


Consider a collection of random variables (not necessarily centralized, nor Nor¬ 
mally distributed), say Yi, Y 2 ,..., and their joint moment-generating function, 
defined by 

M{ti,t 2 ,-) =E[exp(tiYi -^^ 2^2 -f ...)] (4) 

The corresponding joint CUMULANT of Yi, Y 2 , ...Y is defined by 


Ki,2,..J = 


di 

dtidt2---dti 


In [M(ti,t 2 , •■•)] 


tl—*2 


.—t£—0 


(5) 


It is well known (and easy to derive) that = E (Y) , and that all higher- 
order cumulants can be expressed in terms of moments (of the Y variables) 
thus; 

«:i,2,...^= - l)'A‘ji/^j2-Mj. (6) 

jl j2..jieA(l,2,...^) 

where A{1,2, ...t) is the collection of all partitions of the 1,2,...£ indices. A 
PARTITION is a division of a set into an arbitrary number of non-empty and non¬ 
overlapping subsets (these are denoted ji,{ 2 ,...). For example, A{1, 2, 3) consists 
of ji = {1,2,3}, ji,j 2 = {1},{2,3}, ji,j 2 = {2},{1,3}, ji,j 2 = {3},{1,2}, and 
ji,{ 2 , js = {!}, {2}, {3}. 
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There are two points to make about ©: 

• the formula is correct regardless of whether the moments are central or 
simple (from now on, we denote central moments by /ij and simple mo¬ 
ments by fly respectively), 

• using central moments simplifies the RHS substantially - any partition 
containing a single index can be omitted (the corresponding /Tj is zero). 

We will spell out explicitly the first few of these formulas, first using simple 
moments 


1 ^ 1,2 — “ M{1}M{2} 

Ki, 2,3 = ^{1,2,3} “ /^{l,2}/^{3} “ A{1,3}A{2} “ M{2.3}A{1} + ‘^lj-{l}t<'{2}f^{3} 


then, using central moments (note the simplification): 

«1,2 = M{1,2} (7) 

Ki,2,3 = M{1,2,3} 

Ki,2,3,4 = M{1,2,3,4} ~ M{1,2}M{3,4} ~ M{1,3}M{2,4} “ M{1,4}M{2,3} 

To continue, we consider only the special case of having all indices identical (the 
corresponding general formulas get too lengthy): 

Ki,1,1,1,1 = M{i,i,i,i,i} - 10 M{i,i,i}M{i,i} 

«^i, 1,1,1,1,1 = M{i, 1,1,1,1,1} ~ j -I-30/i|]^ ]^| 


Similarly to moments, cumulants are fully symmetric in the permutation of 
indices, e.g. ki_ 2 , 3,4 is the same as k , a , i ,3,2 etc. 

We now proceed to derive explicit formulas for these cumulants when some 
of the Y variables are centralized. Normally distributed (the X’s of the previous 
section - we call them 'singlets ’), and the others are products of two such X’s 
(doublets). This poses a bit of a notational challenge; we will use simple in¬ 
dices for singlets, two indices in parentheses for doublets. For example, Ki, 2 ,( 3 , 4 ) 
indicates a third-order cumulant of three random variables, Xi, X 2 and X 3 X 4 . 

2 Cumulants involving singlets and/or doublets 

It is well known and easy to derive, by combining © and (|5j) , that all cumulants 
involving only singlets are equal to zero, with the exception of 

Hi ,2 = ki,2 
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For doublets, one can derive (by combining ([5]) and (HI), and using a routine, 
‘brute-force’ computation - see the Appendix), that 


K(1,2) 

^(1,2),(3,4) 
^(1,2),(3,4),(5,6) 


fd,2 

^^, 3 ^ 2,4 + ^14X4,3 

^2,3^4,5^1,6 + b2,4V3,5Vi,6 + ^2, 3 ^ 4 , 614 , 5 + b2,4V3_6bd,5 + 

^2,5^3,6^1,4 + V2,6h3,5Vi4 + V 2 ,5^4,6^1,3 + V2,6b4,5Vi_3 


and 


'^(l,2),...(2fc-l,2fc) 


^ ^ ^ 1*2 ^ 3*4 

ie'p2fe 


( 8 ) 


in general, where V 2 k indicates the set of all index pairings (the old V 2 k) which 
do not contain any of the original {1,2} or {3,4} or ... {2k — 1, 2k} pairs. There 
is a simple scheme for building T’ 2 k- 


• Start with the original pairs, e.g. {1, 2}, {3,4}, {5,6}, 

• keeping the first pair fixed, go over all {k — 1)! permutations of the re¬ 
maining pairs, e.g. {1, 2}, {5,6}, {3,4}, 

• keeping the first pair fixed, go over all 2^“^ interchanges of indices in the 
remaining pairs, e.g. {1, 2}, {6, 5}, {3,4}, etc. (four of them in this case), 

• shift each resulting arrangement by one index, e.g. {2,6}, {5,3}, {4,1}. 


One can thus see that the number of terms on the RHS of ([8]) is (fc—1)! x2^“^, 
which equals to 2, 8, 48 when k = 2, 3 and 4 respectively (a fast-growing 
sequence). 


2.1 Mixed cases 

Let us now investigate cumulants with a mixture of singlets and doublets. The 
rules for computing these prove to be quite simple: 

A cumulant involving 

• one singlet (regardless of the number of doublets) equals to zero, 

• more than two singlets (and any number of doublets) is also equal to zero. 

• When a cumulant contains two singlets, it equals to the cumulant in which 
the singlets are replaced by one doublet, i.e. 

^l,2,(3,4),...(2fe-l,2fe) = ^(l,2),(3,4),...(2fc-l,2fc) 
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3 Beyond doublets 

The same approach enables us to develop formulas for cumulants which may 
also involve triplets, quadruplets, etc. These are not needed when dealing with 
parameter estimation related to AR(k) models, but may have a potential appli¬ 
cation elsewhere. 

Thus, for example, odd-order cumulants involving only triplets are all equal 
to zero; for even-orders we get 

«(1,2,3),(4,5,6) = ^{1,2,34,5,6} 

^(1,2,3), (4,5,6), (7,8,9), (10,11,12) = A{1,2,3,4,5,6,7,8,9,10,11,12} 
~A{1,2,3,4,5,6}A{7,8,9,10,11,12} ~ A{1,2,3,7,8,9}A{4,5,6,10,11,12} 
~A{1,2,3,10,11,12}A{4,5,6,7,8,9} 


following 0. Note that the first of these cumulants turns out to be a sum of 
15 terms of the 

bii ,l2 ^3 ,^4 ,^2fc 

type, while the second one has already 9720 terms! 

Conjecture 1 It seems that each cumulant with all indices distinct (regardless 
of how they are grouped) can be always expressed in a form of the RHS of (0!, 
where V 2 k is a specific subset of V 2 k (and 2k is the total number of indices, 
which must be even for all non-zero cumulants). 

Similarly, this is how one evaluates the first few cumulants involving only 
quadruplets: 

^(1,2,3,4) = A{1,2,3,4} 

'^(1,2,3,4),(5,6,7,8) = M{1,2,3,4,5,6,7,8} ~ A{1,2,3,4}M{5,6,7,8} 

«(1,2,3,4),(5,6,7,8),(9,10,11,12) = Ajl,2,3,4,5,6,7,8,9,10,11,12} 
“M{1,2,3,4,5,6,7,8}M{9,10,11,12} ~ A{1,2,3,4,9,10,11,12}A{5,6,7,8} 
“A{5,6,7,8,9,10,11,12}A{1,2,3,4} + 2A{1,2,3,4}M{5,6,7,8}A{9,10,11,12} 

resulting in 3, 96 and 9504 terms, respectively. 

One can then deal with ‘mixed’ cumulants in the same manner. There is of 
course no limit as to their total number; we will give just one example: 

^1,2,(3,4,5),(6,7,8),(9,10,11,12) = A{1,2,3,4,5,6,7,8,9,10,11,12} “ A{1,2}A{3,4,5,6,7,8,9,10,11,12} 
“A{1,3,4,5}A{2,6,7,8,9,10,11,12} “ A{2,3,4,5}A{1,6,7,8,9,10,11,12} ~ M{1,6,7,8}A{2,3,4,5,9,10,11,12} 
“M{2,6,7,8}A{1,3,4,5,9,10,11,12} ~ A{1,2,3,4,5,6,7,8}M{9,10,11,12} ~ M{3,4,5,6,7,8}M{1,2,9,10,11,12} 
+ 2M{1,2}A{3,4,5,6,7,8}M{9,10,11,12} + 2/i{i,3,4,5}M{2,6,7,8}M{9,10,ll,12} 

+ 2M{2,3,4,5}A{1,6,7,8}M{9,10,11,12} 

resulting in a sum of 7848 terms. 
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4 Appendix 

For readers familiar with Mathematica, we now supply a few Mathematica func¬ 
tions to facilitate the computation of cumulants involving centralized, Normally 
distributed random variables and their products. 

4.1 Finding moments 

Computing the mean value (MV) of a product of centralized, Normally dis¬ 
tributed random variables can be achieved by: 

MV[a__] := Module[{aO = Sort[Flatten[{a}]], al, a2, a3}, 
al = Length[aO]; a2 = Table[{aO[[l]], aO[[i]]}, {i, 2, al}]; 
a3 = Table [Delete [aO, {{1}, {i}}], {i, 2, al}]; 

Which[Mod[al, 2] == 1, 0, al == 2, Apply[V, aO], True, 

Apply[Plus, Table[Apply[V, a2[[i]]] MV[a3[[i]]], {i, al - 1}]]]// Expand] 

To use it, we type 

MV[2, 5, 2, 5, 2, 8]/.V[i_, j-]->Vi,j 

where the integers in square brackets represent the corresponding indices, listed 
in any order and allowing for any amount of duplication. 

This returns the following answer: 

3V2,214^512 ,8 + 614^512,8 + 612,212,514,8 


4.2 Computing cumulants 

This time our goal is a bit more ambitious: we want to build a Mathematica 
program for computing cumulants of any number of singlets, doublets, triplets, 
etc. of centralized. Normally distributed random variables. 

We start by constructing a few auxiliary Mathematica functions: 

auxl[b_] := Module[{bl = Union[b], b2, b3}, 
b2 = Map[Position[b, #][[!, 1]] &, bl]; 

b3 = Length[b2]; Table[{b[[b2[[i]]]], Delete[b, b2[[i]]]}, {i, b3}]] 

aux2[a_, B_] := Module[{bl = Map[Join[{a[[l]]}, #] &, 

KSubsets[Drop[a, 1], B[[l]] - 1]]}, 

Flatten[Map[aux3[a, #, B[[2]]] &, bl], Ij] 

aux3[a_, b_, B_] := Map[Join[{b}, #] &, aux4[Complement[a, b], Bj] 

aux4[a_, b_] := Module[{B = auxl[b]}, lf[Length[b] == 1, {{a}}, 

Apply[Join, Table[aux2[a, B[[i]]], {i, Length[B]}]]]] 
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The program to compute a cumulant of any number of products then looks 
like this: 

K[a__] := Module[{aO = {a}, L, q, w}, L = Length[aO]; 
q = Map [Reverse, Partitions[L]]; 

w = Map[aux4[Range[L], #] &, q] /. Table[i -> aO[[i]], {i, L}]; 

Applyplus, Apply[Times, Apply[MV, w, {3}], {2}], {!}] . 
Map[(-l)^®''®*^[^l‘^(Length[^] - 1)! &, q] // Expand] 

To find the fifth-order cumulant of the random variables X 3 , X 1 X 3 , XiX^, 
X 1 X 2 X 3 and XiX 2 X^, one has to type (to simplify the answer, we have assumed 
that the X’s are standardized, i.e. have a mean of zero and the variance equal 
to 1): 


K[3,{l,3},{l,3},{l,2,3},{l,2,3,3}]/.V[L,i_]->l/.V[i_,j.]->Ci,j 


resulting in 

42 -h I 58 CI 2 + 438(71% + 240(7% -h 1784(7i,2(7i,3(72.3 + 1960(7i,2(7%(72.3 
-h802(7%(7% -h 1616(7%C'% -h 240(7% -h 400(7%(7% 

where Cij is the correlation coefficient between Xi and Xj. Note that 


K[{3,l},{3,2,l},3,{2,3,3,l},{l,3}]/.V[L,i_]->l/.V[i_,j.]->Ci,j 


would have returned the same answer. 
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